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Abstract 

We study the crossover from the ordinary to the normal surface universality class in the 
three-dimensional Ising bulk universality class. This crossover is relevant for the behavior 
of films of binary mixtures near the demixing point and a weak adsorption at one or both 
surfaces. We perform Monte Carlo simulations of the improved Blume-Capel model on 
the simple cubic lattice. We consider systems with film geometry, where various boundary 
conditions are applied. We discuss corrections to scaling that are caused by the surfaces 
and their relation with the so called extrapolation length. To this end we analyze the 
behavior of the magnetization profile near the surfaces of films. We obtain an accurate 
estimate of the renormalization group exponent y^^ = 0.7249(6) for the ordinary surface 
universality class. Next we study the thermodynamic Casimir force in the crossover region 
from the ordinary to the normal surface universality class. To this end, we compute the 
Taylor-expansion of the crossover finite size scaling function up to the second order in 
hi around hi = 0, where hi is the external field at one of the surfaces. We check the 
range of applicability of the Taylor-expansion by simulating at finite values of hi. Finally 
we study the approach to the strong adsorption limit /ii — )■ oo. Our results confirm the 
qualitative picture that emerges from exact calculations for stripes of the two-dimensional 
Ising model, [D. B. Abraham and A. Maciolek, Phys. Rev. Lett. 105, 055701 (2010)], 
mean-field calculations and preliminary Monte Carlo simulations of the Ising model on 
the simple cubic lattice, [T. F. Mohry et al, Phys. Rev. E 81, 061117 (2010)]: For certain 
choices of hi and the thickness of the film, the thermodynamic Casimir force changes sign 
as a function of the temperature and for certain choices of the temperature and hi, it also 
changes sign as a function of the thickness of the film. 
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I. INTRODUCTION 



In 1978 Fisher and de Gennes [1[ realized that when thermal fluctuations are 
restricted by a container, a force acts on its walls. Since this effect is analogous 
to the Casimir effect j^, where the restriction of quantum fluctuations induces a 
force, it is called "thermodynamic" Casimir effect. Since thermal fluctuations only 
extend to large scales in the neighborhood of continuous phase transitions it is also 
called "critical" Casimir effect. Recently this force could be detected for various 
experimental systems and quantitative predictions could be obtained from Monte 
Carlo simulations of spin models 0. 

The behavior of the thermodynamic Casimir force can be described by finite size 
scaling (FSS) [^] laws. For the film geometry that we consider here, one gets ^ for 
the thermodynamic Casimir force per area 

Fcasimir — kj^TL^ d(UCi,UC2) 

where Lq is the thickness of the film and t = {T — Tc)/Tc is the reduced temperature 
and Tc the critical temperature. Note that below, analyzing our data, we shall use 
for simplicity the definition t = f3c — /3, where /3 = l/ksT. The amplitude of the 
correlation length ^ is defined by 

^ = ^o,±\try<{l + a±\tr + ct+...) (2) 

where — and + indicate the high and the low temperature phase, respectively. Since 
the correlation length can be determined more accurately in the high temperature 
phase than in the low temperature phase, we take = ^o,+ in eq. ([T]). The power 
law ([2]) is subject to confluent corrections, such as a±|t|'^'^, and non-confluent ones 
such as ct. Critical exponents such as z/ and ratios of amplitudes such as ^o,+/^o,- are 
universal. Also correction exponents such as u and ratios of correction amplitudes 
such as a^/a^ are universal. For the three-dimensional Ising universality class, 
which is considered here uu ~ 0.5. For reviews on critical phenomena and their 
modern theory, i.e., the Renormalization Group (RG) see, e.g., 0-|9j. The universal 
finite size scaling function 6(jjc^^uc2) depends on the universality class of the bulk 
system as well as the surface universality classes U Ci and UC2 of the two surfaces 
of the film. For reviews on surface critical phenomena see e.g. We shall 

give a brief discussion below in section IIIII 

In the past few years there has been great interest in the crossover behaviors of 
the thermodynamic Casimir force. In [l3| the authors have studied the crossover 
from the special surface universality class to the ordinary one by using field theoretic 
methods. They find that for certain choices of the parameters, the thermodynamic 
Casimir force changes sign with a varying thickness of the film. The authors of 



14j | have computed exactly the thermodynamic Casimir force for stripes of the 
two-dimensional Ising model as a function of the external surface fields hi and /i2. 
Also here the authors have found that for certain choices of the fields hi and /;,2, 
the thermodynamic Casimir force does change sign as a function of the tempera- 
ture or the thickness of the film. More recently, the authors of (Tsf have studied 
the crossover from the ordinary to the normal surface universality class, and the 
crossover from the special to the ordinary as well as the normal surface universality 
class using the mean-field approximation. Also in these cases a change of sign of the 
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thermodynamic Casimir force could be observed. Furthermore in isf prehminary 
results 16j of Monte Carlo simulations of the spin-1/2 Ising model on the simple 
cubic lattice for the crossover from the ordinary to the normal surface universality 
class were presented. Following the authors of [l5[ these observations might be of 
technological relevance. They write: "Such a tunability of critical Casimir forces 
towards repulsion might be relevant for micro- and nano-electromechanical systems 
in order to prevent stiction due to the omnipresent attractive quantum mechanical 
Casimir forces j2|, [13]." In recent experiments on colloidal particles immersed in a 
binary mixture of fluids [l8|, the authors have demonstrated that the adsorption 
strength can be varied continuously by a chemical modification of the surfaces. In 
particular the situation of effectively equal adsorption strengths for the two fluids 
can be reached. For sufficiently small ordering interaction at the surface, this cor- 
responds to the ordinary surface universality class. Hence these experiments open 
the way to study the crossover from the ordinary to the normal universality class. 
As discussed in refs. [l9-22| effectively weak adsorption can also be obtained by 
using patterned substrates. 

In the present work we compute scaling functions for the film or plate-plate 
geometry. In order to compare with experiments on the thermodynamic Casimir 
force between colloidal particles and a fiat substrate as studied in ref. [l^ the 
scaling function for the plate-sphere geometry has to be computed. The Derjaguin 
approximation 23| might be used to derive scaling functions for the plate-sphere 
geometry from those for the plate-plate geometry if the radius of the sphere is large 
compared with the distance between the plate and the sphere 24, 2^, as it is indeed 
the case in ref. 18 1. In the recent works j26|, 13] the Derjaguin approximation had 
been used to obtain the scaling functions for the plate-sphere geometry in the strong 
adsorption limit starting from the Monte Carlo estimates of refs. 28, 2^ for the 
film geometry. 

As in ref. ^oj, where we had studied the strong adsorption limit, we shall study 
the crossover by performing Monte Carlo simulations of the improved Blume-Capel 
model on the simple cubic lattice. We shall give the definition of this model in 
section [TTl below. Improved means that corrections to finite size scaling that are 
oc Lq^ vanish. This property is very useful in the study of films, since typically 
the surfaces cause corrections oc 



^0^ 
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and fitting Monte Carlo data, it is 
quite difficult to disentangle corrections that have similar exponents. Motivated 
by the experiments [isj, we shall mainly study films where the external field hi 
at the first surface is finite, while at the other surface the limit /i2 — >■ c>o is taken, 
corresponding to the strong adsorption limit in a binary mixture. For this choice 
of boundary conditions the correlation length of the film divided by its thickness 
remains small at any temperature. In contrast, for /ii = /i2 = the film undergoes a 
second order phase transition in the universality class of the two-dimensional Ising 
model. This implies that in the neighborhood of this transition the correlation 
length of the film divided by its thickness is large. Therefore the Monte Carlo 
study of the crossover from /ii = /i2 = to the limit I/12I — oo would be more 
involved than that performed here. 

In preparation for our study of the thermodynamic Casimir force, we have ac- 
curately determined the surface critical exponent t/h^ of the ordinary surface uni- 
versality class. Furthermore we have estimated the so called extrapolation length 
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for various boundary conditions. The extrapolation length is directly related to 
the corrections to finite size scaling that are caused by the surfaces of the film. 
Our numerical results are mainly based on the analysis of the behavior of the mag- 
netization profile at the bulk critical temperature. Next we have computed the 
thermodynamic Casimir force for the range of inverse temperatures around the 
bulk critical point where, at the level of our numerical accuracy, it is non- vanishing. 



To this end we follow the suggestion of Hucht j3l|]. For alternative methods see 
0, 0, 32, 33[. Note that the stress tensor method of 34] can only be applied for 
periodic or anti-periodic boundary conditions. First we have simulated films with 
a vanishing surface field hi = 0. Based on the data obtained from these simula- 
tions, we have also computed the Taylor-expansion of the thermodynamic Casimir 
force per area in hi up to the second order around hi = 0. We demonstrate that, 
taking into account corrections oc -^^q ^, already for the relatively small thicknesses 
Lq = 8.5, 12.5, and 16.5 the behavior of the thermodynamic Casimir force per area 
as well as its partial derivatives with respect to hi is well described by universal 
FSS functions. Next we have simulated films with various finite values of hi to 
check the range of applicability of the Taylor-expansion and to study the crossover 
beyond this range. Finally we have studied the approach to the strong adsorption 
limit hi — )■ oo. Qualitatively we confirm the picture that emerges from the exact 
solution of the two-dimensional Ising model l3| and the mean- field calculation [l5| . 

The outline of the paper is the following: In section [TT] we define the model and 
the observables that we have studied. In section IIIII we briefiy review the phase 
diagram of a semi-infinite system. Then in section IIVI we discuss the finite size 
scaling behavior of the magnetization profile at the bulk critical point and the 
finite size scaling behavior of the thermodynamic Casimir force. In section |V] we 
discuss how to compute the thermodynamic Casimir force and its partial derivatives 
with respect to the external field hi at the surface. In section I VI I we present the 
results of our Monte Carlo simulations. We performed a series of simulations at 
the bulk critical point, where we focus on the magnetization profile. Next we have 
determined the thermodynamic Casimir force per area in the neighborhood of the 
bulk critical point for various values of the external field hi at the surface. Finally, 
in section IVlII we summarize and conclude. 



II. THE MODEL AND BULK OBSERVABLES 

We study the Blume-Capel model on the simple cubic lattice. It is characterized 
by the reduced Hamiltonian 

H = -P J2 s^Sy + DY^sl-hY^s, (3) 

<xy> X X 

where x = {xo,Xi,X2) denotes a site of the lattice. The components Xq ,Xi and X2 
take integer values. The spin Sx might take the values —1, or 1. In the following 
we shall consider a vanishing external field h = throughout. The parameter D 
controls the density of vacancies Sx = 0. In the limit D ^ —oc the spin- 1/2 
Ising model is recovered. For —oo < D < Dtri the model undergoes a second 
order phase transition in the three-dimensional Ising universality class. For D > 
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Dtri the transition is of first order. Tlie most recent estimate for the tri-critical 
point is Dtri = 2.0313(4) |35|. Numerically, using Monte Carlo simulations it has 
been shown that there is a point {D*,Pc{D*)) on the line of second order phase 
transitions, where the amplitude of leading corrections to scaling vanishes. Our most 
recent estimate is D* = 0.656(20) j36[. In [36[ we have simulated the model at -D = 
0.655 close to (3c on lattices of a linear size up to L = 360. From a standard finite 
size scaling analysis of phenomenological couplings such as the Binder cumulant we 
find 

/3c(0.655) = 0.387721735(25) (4) 

for the inverse of the critical temperature at = 0.655. The amplitude of leading 
corrections to scaling at D = 0.655 is at least by a factor of 30 smaller than for the 
spin- 1/2 Ising model. 

Our recent estimates for bulk critical exponents in the three-dimensional Ising 



universality class are [36 



u = 0.63002(10) , (5) 
?7 = 0.03627(10) , (6) 
w = 0.832(6) . (7) 

In the following we set the scale by using the second moment correlation length 
^2nd in the high temperature phase of the model. On a finite lattice of the linear 
size L in each of the directions it might be defined by 



s2nd 



4sin^ n/L 



where 



(9) 



is the Fourier transform of the correlation function at the lowest non-zero momen- 
tum and 

1 / /V^ 

l3 



x = -{ lE^- (10) 



is the magnetic susceptibility. In [30|, l37|] we find 



6nd,o,+ = 0.2282(2) - 1.8 x (z/ - 0.63002) + 250 x ((3^ - 0.387721735) (11) 

for the amplitude of the second moment correlation length in the high temperature 
phase, where we have used 

t = (3c-(3 (12) 

as definition of the reduced temperature. We shall use this definition of t also in 
the following. The energy density is defined by 

Ebulk = ^ (SxSy) . (13) 

<xy> 
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In the following we shall need the energy density of the bulk system in a neighbor- 
hood of the bulk critical point. To this end, we have performed simulations at 350 
different values of /3 in the range 0.25 < /3 < 0.6 js^. In a small neighborhood of 
/3c, where no direct simulations are available we use 



For a discussion see section IV of 37 . 



A. Film geometry and boundary conditions 

Here we study systems with a film geometry. In the ideal case this means that 
the system has a finite thickness Lq, while in the other two directions the thermo- 
dynamic limit Li, L2 — )■ oo is taken. In our Monte Carlo simulations we shall study 
lattices with Lq <^ Li = L2 = L and apply periodic boundary conditions in the 1 
and 2 directions. 

The reduced Hamiltonian of the Blume-Capel model with film geometry is 
H = -f3j2 s,Sy + Dj24 (15) 

<xy> X 

-Pi ^ s^Sy -132 ^ s^Sy - hi ^ - /i2 

<xy>,xo=yo=i <xy>,xo=yo=Lo x,xo=l x,xo=Lo 

where hi,h2 7^ break the symmetry at the surfaces that are located at xq = 
1 and Xq = Lq, respectively. In our convention < xy > runs over all pairs of 
nearest neighbor sites with fluctuating spins. Note that here the sites (1, xi, X2) and 
(Lq, Xi, X2) are not nearest neighbors as it would be the case for periodic boundary 
conditions. In our study, we set (3i = (32 = throughout. Hence there is no 
enhancement of the coupling at the surface. There is ambiguity, where one puts the 
boundaries and how the thickness of the film is precisely defined. Here we follow 
the convention that Lq gives the number of layers with fluctuating spins. In our 



previous work [30| we have studied the limit of strong adsorption, \hi\, \h2\ — )■ 00. In 
this limit the spins at the boundary are fixed to either —1 or Therefore we had 
put the fixed spins on xq = and xq = Lq + 1 to get Lq layers with fluctuating spins. 
Note that these fixed spins could also be interpreted as external fields hi^2 = ±/3 
acting on the spins at xq = 1 and xq = Lq, respectively. In the following we shall 
denote the type of boundary conditions by (/ii,/i2). In the literature the cases 
/ii = or /i2 = are often called free boundary conditions. To be consistent 
with the literature, we shall denote the strong adsorption limit by + or — in the 
following. In particular the two cases studied in [30'] are denoted by (+, +) = {(3, (3) 
and (+, — ) = (/3, -(3). For the discussion of the behavior of physical quantities near 
the boundary it is useful to define the distance from the boundary. To this end we 
shall assume that the first boundary is located at xq = 1/2 and the second one at 
Xo = Lo + 1/2. Hence the distance from the first boundary is given by z = xq — 1/2 
and the distance from the second one hj z = — xq + + 1/2. 

In order to determine the thermodynamic Casimir force we have measured the 
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energy per area of the film. It is given by 




Since the film is invariant under translations in 1 and 2 directions but not in 
direction, the magnetization depends on xq. Therefore we define the magnetization 
of a slice by 



m(xo) = 




III. PHASE DIAGRAM OF A SEMI-INFINITE SYSTEM 

Here we briefly recall the pha se diagram of a semi-infinite Ising system as it is 



discussed e.g. in the reviews [10Hl2|. For the Blume-Capel model, we expect that 
for D 2 the qual itative features of the phase diagram remain unchanged since 
Dfri = 1.966(2) j38| for the two-dimensional system and Dtri = 2.0313(4) jssf for 
the three-dimensional one. 

In figure [1] we have sketched the phase diagram for a vanishing external field 
h = and a vanishing surface field hi = 0. For P > Pc the spins in the bulk are 
ordered. As a consequence, also the spins at the surface are ordered. This phase 
is denoted by C in figure [H At vanishing bulk coupling /3 = the spins at the 
surface decouple completely from those of the bulk. Hence a two dimensional Ising 
or Blume-Capel model remains that undergoes a phase transition at (3i = /9c,2_d- 
Starting from the point (0, (3c,2d) there is a line of transitions, where the spins at 
the surface order, while those of the bulk remain disordered. This line hits the 
vertical line at /3 = /3c in the so called special or surface-bulk point that we denote 
by SB in figure [H which is a tri-critical point. In figure [H the phase, where both 
the boundary spins and those of the bulk are disordered is denoted by A while the 
one with disordered bulk and ordered surface is denoted by B. The transitions from 
phase A to phase C are so called ordinary transitions, while those from phase B to 
phase C are so called extraordinary transitions. The transitions from phase A to B 
are so called surface transitions. 

For hi the spins at the surface are ordered also for f3i < /3i In the literature 
the transitions from disordered to ordered spins in the bulk for hi ^ are called 



normal transitions. In [39| it has been shown that the normal surface universality 
class is equivalent to the extraordinary surface universality class. 

At the ordinary transition the external field hi at the surface is a relevant per- 
turbation. Hence the RG-exponent i/h^ associated with the surface field is positive. 
In the literature, a number of surface critical exponents have been introduced. In 
the case of the ordinary transition, these can be obtained from and the bulk 
RG-exponents yt = l/i' and yh = {d + 2 — r])/2 hj using scaling relations. In the 
following we need 

^l = mn , (18) 

f3i = iy{d-l-yf,,) , (19) 
ji = u{2-d/2-r]/2 + yh,) . (20) 
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FIG. 1. Sketch of the phase diagram of the semi-infinite system. On the x-axis we plot 
the couphng (3 of the bulk and on the y-axis the excess coupling /3i of the surface. A 
detailed discussion is given in the text. 



For the definitions and a complete list of these exponents see the reviews |10l-ll2 . 
The numerical values of surface critical exponents for the three-dimensional Ising 
universality class have been computed by various theoretical methods. Mean field 
theory predicts = 1/2. The authors of |40[| q uote = 0.7363 as result of their 
real space RG method and the authors of [41[ quote 71 = 0.78(2) as result of a 
series expansion, which corresponds to = 0.72(3). The e-expansion gives [42 



1 31 

-e H ( 

6 321 



(21) 



Naively inserting e = 1 one gets yt^ = 0.666... and yh^ = 0.762... at 0(e) and O(e^), 
respectively. Using a massive field theory approach the authors of [43] obtain Ai = 
0.45 from the [1/1] Fade approximant of their two-loop result, which corresponds to 
yhi = 0.714. Comparing the different Fade approximants that are given in table 9 



of [43|] one might conclude that the uncertainty of the estimate of yhi is about 0.02. 
In table [T] we have summarized Monte Carlo results for surface critical exponents. 
Most of the authors quote an estimate for /3i and some in addition for 71. In those 
cases in which the authors did not quote a result for yh^ we have converted the 
value given for (3i using the scaling relation ( IT9|) and u = 0.63002(10). 
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TABLE I. Monte Carlo results for surface critical exponents for the ordinary phase tran- 
sition in the three-dimensional Ising universality class. The authors of [4^^ quote no final 
result for 71. Here we give the average of the three results given in table II of j^j . In 
case the authors do not quote an estimate for y^^, we have computed it from j3i and the 
scaling relation i li9j) . These cases are marked hy * . 

Ref. /?! 71 Vhi 

[44] 0.78(2) 0.762(32)* 

[45] 0.79(2) 0.79(10) 0.746(32)* 

[46] 0.78(2) 0.78(6) 0.762(32)* 

[47] 0.740(15) 

[48] 0.807(4) 0.760(4) 0.719(6)* 

[49] 0.80(1) 0.78(5) 0.730(16)* 

[50] 0.737(5) 

[51] 0.796(1) 0.7374(15) 

[52] 0.795(6) 0.738(10)* 

here 0.7249(6) 



For comparison we also anticipate our result for yh^ that we obtain in section 
IVIAI below. Except for 48| the estimates for are larger than ours. In particular, 
note that the difference between our result and that of 
large as the combined error. 



51 is about six times as 



IV. FINITE SIZE SCALING APPLIED TO FILMS 

In this section we shall discuss the finite size scaling behavior of the magnetiza- 
tion profile at the bulk critical point and thermodynamic Casimir force for arbitrary 
temperature. The starting point of our considerations is the reduced excess free en- 
ergy per area of the film 

fexiLf), t, hi) = ffiim{Lo, t, hi) — Lofbuikit) (22) 

where ffUmiLo, t, hi) is the reduced free energy of the film per area and fbuikit) the 
reduced bulk free energy density. There is no dependence on h2, since we consider 
the limit /12 — ^ 00. The singular part of the reduced excess free energy per area has 



the finite size scaling behavior |1CH12 



fecc,s{Lo,t,hi) = L^''^'g{t[Lo/^or,hi[Lo/le„,or''^) (23) 

where we have ignored corrections to scaling at the moment and = 3 is the 
dimension of the bulk system. We shall define the amplitude lex,nor,o of the normal 
extrapolation length lex,nor below, eq. fl30|) . Note that the bulk contributions to 
the non-singular part of the free energy cancel in eq. fl22|) . However, there remain 
contributions from the two surfaces. 
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A. The magnetization profile at the bulk critical point 



In terms of the reduced free energy per area the magnetization at xq = 1 is given 



by 



mi 



dfex{Lo,t, hi) 
dhi 

L2 Z dhi 




^{l,Xl,X2) 



(24) 



Xl,X2 



In section lVI Al we shall determine the value of the RG-exponent yti from the scaling 
of nil with the thickness Lq at /ii = and (3 = 13^. Taking the partial derivative of 
eq. (125]) with respect to hi we get 



mi 



dfe. 



dgit[Lo/^of\hi[Lo/lex,nor,o\''''^) 



dhi 

T-d+l 
^0 



t=hi=0 



dhi 



t=hi=0 



t=hi=0 



ex, 



nor,0 



(25) 



where denotes the partial derivative of g with respect to = /ii[^vo//ea-,nor,o]^''i- 
Note that the non-singular contribution to fex from the first surface does not feel 
the breaking of the symmetry by the second surface. Therefore it is an even function 
of hi and does not contribute to the partial derivative with respect to hi. 



The extrapolation length l^x can be defined by the behavior |10Hl2 

m{xo) = cLq'^'" ipixo/Lo) 



(26) 



of the magnetization profile at the critical point of the bulk system. Note that from 
scaling relations it follows that P/u = [1 + '7)/2, where r] = 0.03627(10) for the 
three-dimensional Ising universality class jssj. 

In the neighborhood of the surface with spins fixed to = 1, one expects that 
for z Lq, where z = Lq — xq + 1/2, the magnetization profile does not depend on 
Lq. Therefore ^(a^o/^o) = {z/Lq)-^/" and hence [iMl 



m{XQ) 



cz 



-13/u 



(27) 



Also at the free boundary we expect that for z ^ Lq, where now z = xq — 1/2, the 
functional form of the magnetization profile does not depend on Lq. As we have 
seen above, for a fixed value of z, the magnetization behaves as mi oc Lq '^'^^'^^'^^ . 



Therefore 10-12 



m[Xo) 



a z 



-P/u+d-l-yh^ 



a z 



(28) 



Since —fi/v < 0, the scaling function of the magnetization profile diverges as 
z/Lq ^ at the boundary with fixed spins. On the other hand since {(3i — P)/v > 
the scaling function of the magnetization vanishes as z / Lq ^ at the free boundary. 

Based on this observation one might define for finite thicknesses Lq an effective 
distance from the boundary 



Zeff = Z + le 



(29) 
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such that the magnetization profile at Zejj = vanishes for hi = or diverges in the 
case of symmetry breaking boundary conditions. The concept of the extrapolation 
length has been worked out explicitly for the ordinary transition in the framework 
of mean-field theory [13] • Also in the Monte Carlo study of the magnetization 
profile of a semi-infinite system in the extraordinary surface universality class an 
extrapolation length had been introduced j53|. The extrapolation length is related 
with corrections oc Lq^ discussed in the framework of field-theory in [HJ]. In the 
following we shall distinguish between the extrapolation length lex,ord [ord for ordi- 
nary surface transition) and lex,nor {nor for normal surface transition) in the case of 
symmetry breaking boundary conditions. The extrapolation length depends on the 
precise definition of z. Physically, the extrapolation length depends on the details 
of the microscopic model, in particular on the details of the fields and interactions 
at the surface. In section IVIBI we shall study the behavior of the extrapolation 



length as a function of the field hi at the boundary. One expects [55| 

^ex,noir{hl) ^ex,nor,0^l ("^0) 



which defines the amplitude lex,nor,o that we have already used above in eq. fl23|) . 
Capehart and Fisher [5^ have argued that the arbitrariness in the definition 



of the thickness of the film leads to corrections oc Lq^. These corrections can be 
eliminated by replacing Lq in finite size scaling laws such as eq. ( 12^ by an effective 
thickness 

Lo,eff = Lo + Ls (31) 

of the film. Assuming that the corrections due to a surface are caused by a unique 
irrelevant surface scaling field, the constant Lg should be given by 

Ls = hx,l + hx,2 (32) 

where lex,i and lex,2 are the extrapolation lengths at the two surfaces of the film. 
In section II A 4 of ref. 13] a similar discussion of the extrapolation length had 



been presented. For a discussion of the effective thickness and further references 



see section IV of ref. 30 



B. Crossover scaling function of the thermodynamic Casimir force 

In terms of the reduced excess free energy per area the thermodynamic Casimir 
force per area is given by [5| 

-'^Casimir r%j- ■ x'^'^J 



Using the finite size scaling law we arrive at 



dLo 



+ L^'^yhlhilLo/le^^norfiY'"' ghiiXt, XhJ (34) 



where Xt = t[Lo/^o]^* and = hilL^/ lex,nor,o]^''^ ■ The partial derivatives of g with 
respect to Xt and x^.^ are denoted by gt and gh^ , respectively. Note that the analytic 
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part of fex is due to the surfaces and does not depend on Lq and therefore does not 
contribute to the thermodynamic Casimir force. It follows that the thermodynamic 
Casimir force per area follows the finite size scaling law [l5[ 

Fcasimir = keT L^''Q{t[L^ / hilL^/lex^norflY"') (35) 

where 

0(xt,XhJ = {d-l)g{xt, Xh^)-ytt[LQ/ ^q]^' gt{xt, Xh)-yhihi[Lo/lex,nor,of''^ gh^Xt, xh^) ■ 

(36) 

Taking the n^^ derivative of the thermodynamic Casimir force with respect to hi 
we get 

d Fcasimir , rpj-d\T n inVhi 0(^t;^hi) 

-— = kBTL^ [Lo/lex^norfl] ^-^^ • (37) 

1 h\ 



V. COMPUTING THE THERMODYNAMIC CASIMIR FORCE AND 
DERIVATIVES WITH RESPECT TO THE EXTERNAL FIELD AT THE 
SURFACE 

On the lattice, we approximate the derivative of the reduced excess free energy 
per area with respect to the thickness Lq of the film by a finite difference: 

df 

^ A/e.(Lo) = fe.{L, + 1/2) - /e.(Lo - 1/2) (38) 



where Lq + 1/2 and Lq — 1/2 are positive integers. As suggested by Hucht [311 we 
compute this difference of free energies as the integral of the difference of corre- 
sponding internal energies: 

Af,xiLo,P) = - d/3AEe.(Lo,/3) (39) 

where 

AE,,{Lo) = E{Lo + 1/2)- E{Lo- 1/2)- Ebuik • (40) 

In practice the integral (139|) is computed by using the trapezoidal rule. Our previous 
experience [30| shows that AEex{Lo) has to be computed for about 100 values of /3 
to obtain Afex{Lo,P) in the whole range of /3 we are interested in at the level of 
accuracy we are aiming at. 

In this work we compute the Taylor-expansion of the thermodynamic Casimir 
force with respect to the boundary field hi around /ii = up to the second order. 
To this end we compute the first and second derivative of Af^xiLo) with respect to 
hi. The n*'^ derivatives can be written as 

a"A/e.(Lo,/3,/ii) ^ _ ~ d^AE,xiLoJ,hi) 

dh-i i^/^ dh-i ^ ' 

where 

d^AE,x{Lo,P,hi) _ d^E{Lo + l/2,P,hi) d^Ejlp - 1/2, P, hi) 
dhl dhl dhl 
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(42) 



Note that there is no bulk contribution, since the internal energy of the bulk does not 
depend on hi. In the Monte Carlo simulation, the first derivative can be computed 
as 

OEjLo.S.h,) 
dhi 

where 

E 



{EMi)-{E){Mi) (43) 



L2 

<xy> 

and 

Ml = 5^ . (45) 

X\,X2 

The second derivative is given by 
a2^(Lo,/3,/i 



dhl 



^ = {EMl) - 2{EMi){Mi) - {E){Ml) + 2{E){Mif . (46) 



Higher derivatives could be computed in a similar way. However it turns out that 
the relative statistical error of the second derivative is much larger than that of the 
first one. Therefore it seems useless to implement and measure higher derivatives. 



VI. MONTE CARLO SIMULATION 

First we have simulated films with (0, +) boundary conditions at the bulk critical 
point for thicknesses up to Lq — 64. Analyzing the data obtained from these 
simulations, we have determined the value of for these boundary conditions and 
have obtained an accurate result for the RG-exponent yhi- Next we have simulated 
lattices of the size Lq = L = 512 with (+,0) and (/ii,0) boundary conditions 
with hi — 0.2, 0.1, 0.05 and 0.02 at the bulk critical point. Prom the behavior of 
the magnetization profile in the neighborhood of the surfaces we have determined 
the extrapolation length lex,ord for free boundary conditions and the extrapolation 
length lex,nor ^s a function of hi. Next we have studied {hi, — ) boundary conditions 
for hi = 0.2, 0.18, 0.16, 0.15, 0.14, 0.13, 0.12, 0.11, 0.1, 0.09, 0.08, 0.07, 0.06 and 
0.05 also at the bulk critical point. Prom the zero of the magnetization profile, we 
read off the difference lex,norihi) — lex,nori—) of extrapolation lengths. Note that 
lex,nor{-) = lex,nor{+) duc to Symmetry. 

Next we have studied the thermodynamic Casimir force per area in the neigh- 
borhood of the bulk critical point. To this end, we have simulated films of the 
thicknesses Lq = 8, 9, 12, 13, 16 and 17 for about 100 values of /3 each. Using 
the data obtained from these simulations we have computed the finite size scaling 
function of the thermodynamic Casimir force per area for (0, +) boundary condi- 
tions. Purthermore we have computed the Taylor-expansion of the thermodynamic 
Casimir force per area for (/ii,+) boundary conditions to second order around 
hi = 0. We have simulated Lq = 8 and 9 at hi = 0.03, 0.06, 0.1 and 0.2 to check for 
how large values of hi and hence of Xhi the Taylor-expansion accurately describes 
the finite size scaling function Q{xt,Xhi). Finally we have studied the approach to 
the strong adsorption limit as Xhi — > oo. 
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As in our previous work [30[ we have simulated the Blume-Capel model by using a 
hybrid js^l of local heat-bath updates and cluster updates [58, Hj. Since the cluster 



updates only change the sign of the spins, additional local updates are needed to 
ensure ergodicity of the compound algorithm. In one cycle of our algorithm we 
sweep twice through the lattice using the local heat bath algorithm followed by one 
or more cluster-updates. In one sweep we run through the lattice in typewriter 
fashion, performing heat bath updates site by site. We have always performed a 
cluster-update, in which all spins are flipped that are not frozen to the boundary. 
For a detailed discussion see section V A of ref. js^. Note that here, in contrast 
to ref. jsoj, we have applied this type of cluster- up date also to systems with (+, — ) 
boundary conditions. To this end we had to adapt the implementation of the cluster 
search; we had to allow for the possibility that two spins in the cluster frozen to 
the boundary might have different signs. Furthermore, we have generalized the 
cluster-update to the case of a finite external field hi at the surface. A spin at the 
boundary freezes to the external field with the probability pf = 1 — Pd, where 

Prf = min[l,exp(-2/;,iS^)] . (47) 

In the case of large systems, discussed in sections I VI B 1 1 IVI B 21 below, we per- 
formed in addition single cluster updates jsoj. In all our simulations we have used 
the SIMD-oriented Fast Mersenne Twister algorithm ^] as pseudo-random number 
generator. 



A. Simulations at the bulk critical point 

First we have simulated films with (0, -f-) boundary conditions at our estimate 
of the bulk critical point /3c = 0.387721735 [39]. Since the fixed spins at the second 
surface act effectively as an external field for the effectively two-dimensional system, 
the correlation length of the film stays finite at any value of p. This means that for 
a given thickness Lq, finite L effects decay oc exp{—L/C,fiim) for sufficiently large 
values of L. Hence we can chose L such that finite size effects are much smaller 
than the statistical errors and therefore can be ignored in the analysis of our Monte 
Carlo data. In order to check which values of L are needed to this end, we have 
performed simulations for the thickness Lq = 6, using L = 6, 7, 8, 9, 10, 11, 12, 13, 
14, 15, 16, 18, 20, 24, 32 and 48. For each of these lattice sizes we have performed 
10^ or more update cycles. As a check we have simulated films with the thickness 
Lq = 12 and L = 24, 24, 28, 32, 40, 48, 56, 64 and 96, where we performed 10^ 
update cycles throughout. We have studied the behavior of the second moment 
correlation length, the magnetization at the surface mi, and the energy per area of 
the film E and its first and second derivative with respect to hi. 

Here we use the same definition of the second moment correlation length as in ref. 



30| . See in particular section HI C of [30|. The disadvantage of this definition of the 



second moment correlation length is that as soon as more than one eigenstate of the 
transfer matrix contributes to the correlation function, corrections to the L — )■ oo 
limit only decay oc L~^. In figure [2] we have plotted the second moment correlation 
length obtained with the pairs of wave vectors ( (0, 0), (0, 1) ) and ( (1, 0), (1, 1) ) as 
a function of L~^. While the estimate obtained by using the pair of wave vectors 
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FIG. 2. We plot (,2nd as a function of L~'^, where L is the hnear extension in the transversal 
directions, for films of the thickness Lq = 6. We have computed the second moment corre- 
lation length by using the pairs of wave vectors ( (0, 0), (1,0) ) (circles) and ( (1, 0), (1, 1) ) 
(triangles) . 



( (1, 0), (1, 1) ) is monotonically increasing with increasing L, the estimate obtained 
by using the pair ( (0,0), (0, 1) ) displays a minimum close to L = 12. The value 
at this minimum is about 0.993 times the asymptotic value. Fitting the results 
obtained for L = 24, 32 and 48 with the ansatz $,2nd{L) = ^2nd + aL^^ we get 
i2nd = 1.6988(6) and 1.6990(5) for the choices ( (0, 0), (1, 0) ) and ( (1, 0), (1, 1) ), 
respectively. 

Next we have analyzed the energy per area, its first and second derivative with 
respect to hi and the magnetization mi at the surface. These quantities should 
converge with exponentially small corrections as L — oo. We have fitted these 
quantities with the ansatz A{L) = A{oo) + ca exp(— L/^jj/m), where we have taken 
our result for the second moment correlation length C,2nd = 1-70, which should not 
be much smaller than the exponential correlation length that is actually needed 
here. Fitting all data with L > 16 we find for the magnetization at the boundary 
xVDOF = 4.55/4, mi(oo) = 0.1250175(6) and = -0.237(25). This means that 
for L ^ lie, film the deviation from the limit L — )■ oo has about the same size as the 
statistical error that we have reached here for Lq = 6. Note that below, for larger 
thicknesses the number of measurements is more than a factor of ten smaller than 
for Lq = 6. Analyzing the energy per area and its first and second derivative with 
respect to hi we find that for L ^ IOC, fum the deviation from the limit L — )■ cxd has 
about the same size as the statistical error. Analyzing our results for the thickness 
Lq = 12 we find consistently that for mi and the energy per area and its first 
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and second derivative with respect to hi, the deviation from the hmit L — )■ oo has 
about the same size as the statistical error for L ^ 10C,fUm- As we shall see below, 
^fiim ~ 0.225(Lo + 1.43). Therefore, for L = 4Lo, which we have used below, the 
deviation from the limit L — > oo should be clearly smaller than the statistical error 
and can hence be ignored. 

Next we have simulated films for a large number of thicknesses up to Lq = 64 at 
f3 = 0.387721735, using L = ALq throughout. For the thicknesses Lq = 6, 7, 8, 9, 
10, 11, 12, 13, 14, 15, 16, 18, 20, 22, 24 we performed 10^ update cycles throughout, 
and 7.6 x 10^ 10^ 8.7 x 10^ 6.5 x 10^ 4.3 x 10^ 2.6 x 10^ and 2.5 x 10^ update 
cycles for Lq = 28, 32, 36, 40, 48, 56, and 64, respectively. These simulations took 
about 18 months of CPU time on a single core of a Quad-Core AMD Opteron(tm) 
Processor 2378 running at 2.4 GHz. 

We have fitted the data of the magnetization mi at the surface with free boundary 
conditions with the ansatz 

mi = 6 (Lo + (48) 

where b and are the parameters of the fit and 

mi = 6 (Lo + L,)2-J"^i X [1 + c (Lo + L,)-^] (49) 

where now c is an additional parameter, to obtain some control on sub-leading 
corrections. We have taken into account all data obtained for the thicknesses Lq > 
Lo^min- Fitting our data with the ansatz f H5]) we get acceptable fits already for 
Lo,min = 10: h = 1.6131(13), = 1.4289(33), yn^ = 0.72493(20), and xVDOF 
= 13.2/15. Fitting with ansatz (149|) we get for Lo^rnin = 6 the results b = 1.6109(22), 
Ls = 1.4166(86), Vh, = 0.72520(31), c = -0.09(4), and xVDOF = 18.1/18. 
We arrive at the final estimates 

6=1.613(4) (50) 
Ls = 1.43(2) (51) 
Vh, = 0.7249(6) (52) 

where the central result is taken from the fit with the ansatz fHHj) and Lo,mm = 10. 
The error bar is chosen such that also the result for the fit with sub-leading correc- 
tions fH9|) is covered. We have also estimated the error induced by the uncertainty of 
our estimate of the inverse bulk critical temperature Pc- To this end, we have first de- 
termined the derivative of mi with respect to (3 for Lq = 8, 9, 12, 13, 16 and 17, where 
we performed simulations for many values of /3. We have extrapolated these results 
to other values of Lq assuming dnii/dp oc [Lq + L^)^"^**!^^*. Using this we have 
computed mi at /3 = /3c -|- error = 0.38772176 and have redone the fits performed 
above. We find that the deviations of the results for /3 = 0.38772176 from those for 
P = 0.387721735 are much smaller than the errors quoted in eqs. ( !50|51|52p . 

Next we have analyzed the second moment correlation length obtained by using 
the pair ( (1, 0), (1, 1) ) of wave vectors. Following the discussion above, finite L 
effects might be still at the level of 1% for our choice L = ALq. Since this effect is 
essentially the same for all thicknesses, it mainly effects the parameter c in the two 
equations below. First we have fitted our data with the ansatz 

6nd = c (Lo + Ls) (53) 
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where c and Ls are the parameters of the fit. Using Lo^min = 8 we obtain c = 
0.22435(9), Ls = 1.487(6) and xVDOF = 19.2/18. Fitting instead with the ansatz 

6nd = c (Lo + Ls) X [1 + 6 (Lo + Ls)-^] (54) 

we get for Lo,mm = 6 the results c = 0.22476(16), = 1.422(20), b = 0.48(12) and 
X^/DOF = 15.0/19. We conclude that the estimate for Lg obtained from the finite 
size scaling behavior of ^2nd is consistent with but less precise than that obtained 
from the finite size scaling behavior of mi. 

Finally we have fitted the energy per area with the ansatz 

E = L^Ens + Ens,s + c (Lo + L,) "2+1/- (55) 

where we have used Ens = 0.602111(1) [37] and u = 0.63002(10) as input. Starting 
from our smallest thicknesses we get acceptable fits: For -Lo,mm = 6 we obtain 
c = -3.5916(7), Ls = 1.4136(21), Ens,s = 3.0644(2) and xVDOF = 18.8/19. As 
check we have also fitted with the ansatz 

E = LoEns + Ens,s + c (Lq + L,)-2+i/- x [1 + 6 (Lq + Ls^^] (56) 

where we have included sub-leading corrections. For Lo^rnin = 6 we get c = 
-3.5930(15), Ls = 1.423(12), Ens,s = 3.0646(2), b = 0.017(23) and xVDOF 
= 18.4/18. 

We have also redone the fits for shifted values of Ens and u. Since we have seen 
above in the case of mi that the uncertainty of Pc is negligible, we have skipped this 
check here. Taking all these results into account we arrive at the final estimates 

Ens,s = 3.064(1) (57) 
Ls = 1.42(2) (58) 
c = 3.592(3) . (59) 

In particular we notice that the estimate of Ls is fully consistent with that obtained 
above from the analysis of the magnetization mi at the boundary. In the following 
we shall use Ls = 1.43(2) as obtained from the analysis of the magnetization mi at 
the free boundary. 



B. The extrapolation length 

First we have simulated lattices with {hi,0) boundary conditions of the size 
Lo = L = 512 at /3 = 0.387721735 using hi = (3^, 0.2, 0.1, 0.05 and 0.02. For this 
geometry one expects strong finite L effects. However these should not alter the 
behavior in the neighborhood of the boundary that we study here. In all cases we 
have performed 2.6 x 10^ update cycles. In total these simulations took about 7 
months of CPU time on a single core of a Quad-Core AMD Opteron(tm) Processor 
2378 running at 2.4 GHz. 
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1. Behavior of the magnetization at the free boundary 



Following the discussion of section IIV Al we have determined the extrapolation 
length lex,ord by fitting our data for the magnetization profile with the ansatz 

m{z) = c{z + h,,ordf^^-^^''' (60) 

where z gives the distance from the boundary as defined in section III At a few lines 
above eq. (fT6!) . To this end, we have computed the ratios 

r{z) = m{z + I /2)/m{z- 1/2) (61) 



to eliminate the constant in eq. fl60|) . It turns out that cross-correlations of these 
ratios are relatively small. Therefore, for simplicity, we have fitted our data for these 
ratios taking only their statistical error into account, ignoring cross-correlations. 
The statistical errors of the fit-parameters were computed by using a Jackknife 
procedure on top of the whole analysis, providing us with correct statistical errors 
for the results. 

First we have fitted our data with the ansatz 

r(^)=( 't^"'!!/9 )''^"''''^ (62) 

where the free parameters of the fit are the extrapolation length lex,ord and the 
exponent (/3i — /3)/z/. We have performed a large number of fits with various choices 
of the range Zmin < < Zmax of distances from the boundary that are taken into 
account, for all values of hi that we have simulated. The results for different hi 
are consistent among each other. In figure [3] we show our results for the exponent 
{(3i — (3)/u for the choice z^ax = ^Zmin as a function of Zmim where we have averaged 
over all values of hi that we have simulated. The error that we give is purely 
statistical. For comparison we plot the estimate of (/3i — /3)/z/ = 2 — yh^ — {l + ri)/2 = 
0.7570(7) obtained by using our estimate of eq. ( 15^ . and rj = 0.03627(10) js^. 
We find that for Zmin = 5 up to 30 the estimates obtained from the behavior of the 
magnetization profile in the neighborhood of the surface are consistent with but 
less precise than the one using the estimate of i/hi obtained in the previous section. 
Therefore, in order to determine our final result for the extrapolation length lex,ord, 
we have fixed (/3i — /3)/z/ = 0.7570(7). Fitting the data for r{z) averaged over all 
values of hi that we have simulated in the range 5 < z < 30 we arrive at 

lex,ord = 0.48(1) (63) 

where the error is dominated by the uncertainty of {/3i — (5)/v. 



2. Normal extrapolation length as a function of hi : part 1 

Following the discussion of section IIV Al the magnetization in the neighborhood 
of the surface behaves as 

m{z, hi) CX (Z + lex,nor{hi))~^''' (64) 
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FIG. 3. We plot the estimate of (/3i — f3)/v obtained by fitting with ansatz (I62p as a 
function of Zmin (filled circles). In these fits distances Zmin ^ ■z < 'izmin from the boundary 
are taken into account. We have averaged the results over all values of hi that we have 
simulated. These results are compared with (/3i — f])/!^ = 2 — yi — [1 + r])/2 = 0.7570(7) 
obtained by using yh^ = 0.7249(6), see the previous section, and r] = 0.03627(10) |36l ]. 
where the central value is depicted by the solid line and the error-bars are indicated by 
the dashed lines. 



where z gives the distance from the boundary. Also in the case of symmetry break- 
ing boundary conditions, we have computed ratios (1^ of the magnetization of 
neighboring shces. These behave as 

^1^^^ I 'ex, nor -■-/ ^ \ (65) 

\^ ~l~ ^ex,nor 1/^/ 

Here we have solved eq. (!65|) with respect to lex,nor for a single value of z, where we 
have used P/p = 0.518135. For hi = f3c we find for z ~ 15 only a small dependence 
of the result on z. We read off lex,nor = 0.96(2). In a similar way we have determined 
the extrapolation length for the other values of hi. Our results are summarized in 
table HI 

In ref. j36| we had determined Lg = 1.9(1) for (+, +) and (+, — ) boundary 
conditions analyzing films of thicknesses up to Lq = 32 at the critical point of the 
bulk system. Now we have added for (+, +) boundary conditions simulations for 
the thicknesses Lq = 48, 64 and 96. This allows us to improve the accuracy of 
our estimate. Now we find Lg = 1.90(5). This result is in perfect agreement with 
Lg = 2lex,nor{(3c) = 1-92(4) obtained here. 

For (0,+) boundary conditions we find Lg = lex,ord + lex,nor{Pc) = 0.48(1) + 
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TABLE II. The extrapolation length lex,nor is obtained for various values of hi by ana- 
lyzing the behavior of the magnetization profile near the surface. For a discussion see the 
text. 



"-1 '•ex,nor 

Pc 0.96(2) 
0.2 2.25(3) 
0.1 5.56(4) 
0.05 14.0(2) 
0.02 PS 45 



0.96(2) = 1.44(3), which is in perfect agreement with the result given in eq. f lSTj) 
above. 

3. Normal extrapolation length as a function of hi : part 2 

Here we have simulated systems with hih2 < 0, where h2 = —f3c, corresponding 
to fixed spins Sx = — 1 at = Lq + 1, and various values of hi. For such a 
choice of boundary conditions the magnetization profile takes positive values in 
the neighborhood of the first surface and negative ones in the neighborhood of the 
second surface. Therefore, in between the magnetization profile vanishes at Xo^zero, 
which depends on hi and /i2- The distance of this zero from the first boundary is 
given by Xo,^ero — 1/2 and from the second boundary by Lq — a;o,^ero + 1/2. Our basic 
assumption is that the zero of the magnetization indicates the physical middle of 
the system. Hence the distances of the zero from the effective positions of the first 
and the second boundary should be the same: 

^O,zero ~ 1/2 + lex,nor{hi) = Lq — XQ^zero + 1/2 + lex,nor{h2) (66) 

and hence 

^lex,nor{hi, /i2) = hx,nor{hl) — hx,nor{h2) = Lq + 1 — 2XQ^zero ■ (67) 

In order to define the zero of the magnetization we have linearly interpolated the 
magnetization profile. Throughout we simulate lattices with L = 4Lo. First we have 
simulated at hi = 0.2, 0.1 and 0.05, using a large number of thicknesses Lq. Our re- 
sults are summarized in table IIIII Apparently, Alex,nor converges with an increasing 
thickness Lq. Numerically, corrections to the Lq — )■ oo limit are compatible with an 
exponential decay. However we can not strictly exclude power-like corrections. For 
hi = 0.2 our results for Lq > 20 are compatible within the statistical error. In the 
case of /ii = 0.1 the results for Lq = 40 and 48 are compatible. The one for Lq = 64 
is larger by about twice the combined statistical error than that for Lq = 48. For 
hi = 0.05 the results for Lq = 120 and 160 are compatible. It is natural to assume 
that the thickness Lq needed to obtain Alex,nor with a given relative error is propor- 
tional to the extrapolation length lex,nor{hi). Using lex,nor{(3c) = 0.96(2) obtained 
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in the section above, we conclude that for Lq ~ I0lex,nor the deviation of Alex,nor 
from its Lq — )■ oo hmit is less than the statistical error that we have reached here. 
Next we have simulated at h = 0.18, 0.16, 0.15, 0.14, 0.13, 0.12, 0.11, 0.09, 0.08, 
0.07 and 0.06 for a single thickness Lq each. The thicknesses Lq and the estimates 
for AZgx.nor are given in table UTTl Throughout Lq > 10lex,nor holds. For each of the 
simulations given in table UTTl we performed about 10^ update cycles. In total these 
simulations took about 8 months of CPU time on a single core of a Quad-Core AMD 
Opteron(tm) Processor 2378 running at 2.4 GHz. Note that the results obtained 
here, are consistent with those of the previous section. Taking the numbers from 
table [HI we get A/ex,nor = 1.29(5), 4.60(6) and 13.04(22) for hi = 0.2, 0.1 and 0.05, 
which is perfectly consistent with the results of the present section, given in table 

m 

From eq. (|30|1 follows that 

Alex,nor = Iq + lex,nor,o\hl\'^^^''^ (68) 

where naively Iq = lex,nor{Pc)- However, since lex,nor depends on the precise definition 
of the thickness of the lattice, we keep the offset Iq as a free parameter here. 

It turns out that for the range of hi that we have simulated here, analytic 
corrections have to be taken into account. Therefore we have fitted our results for 
the difference of the extrapolation length with the ansatz 

Alex,nor = k + lex,nor,o\hl + fl/lil"^/^"! (69) 

where the amplitude lex,nor,o , the offset Iq and the correction amplitude a are the 
parameters of the fit. Note that there should be no term oc hf since lex,nor{hi) = 
hx,nor{~hi). Wc sct i/h-^ = 0.7249(6) as obtained in section IVI A[ In addition we 
have fitted with 

Alex,nor = lo + lex,nor,o\hl + ah\ + | "^/^''^i (70) 

to check for systematic errors due to the truncation of the Wegner expansion. Al- 
ternatively we have also fitted with 

A/ex,nor = ^0 + lex,nor,o\hl\~^'^''^ X (1 + S/l?) (71) 

and 

A/ex,nor = ^0 + 4x,nor,0 1 /^l T^^^"! X (1 + o/l^ + hh{) . (72) 

Fitting with the ansaetze fl69|7ip we get acceptable values of x^/DOF starting from 
hi^max = 0.2, i.e. taking all data into account. Discarding data with large hi 
the result for lex,norfi is slightly decreasing and also x^/DOF is further decreasing. 
E.g. fitting with ansatz ( E^ and taking hi ^ax = 0.14 we get lex,norfi = 0.2133(9), 
Iq = 0.04(14), a = 6.9(2.2) and xV^OF = 1.72/7. Taking into account the variation 
of the results over various ansaetze that we have used and the uncertainty of i/hi 
we arrive at 

lex,nor,0 = 0.213(3) (73) 

which we shall use in the following. 
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TABLE III. The difference of tlie extrapolation lengtlis Alex,nor{hi, = lex,nor{hi) — 
hx,nor{h2), where h2 = —Pc as a function of hi. For a discussion see the text. 



hi 


Lo 




0.2 


10 


1.2642(21) 


0.2 


12 


1.2748(25) 


0.2 


14 


1.2829(27) 


0.2 


16 


1.2889(32) 


0.2 


18 


1.2852(35) 


0.2 


20 


1.2955(38) 


0.2 


22 


1.2989(42) 


0.2 


24 


1.2938(44) 


0.2 


28 


1.2932(54) 


0.2 


32 


1.2941(62) 


0.18 28 


1.6182(47) 


0.1(3 32 


2.043-1(61) 


0.15 36 


2.313(7) 


0.14 42 


2.606(9) 


0.13 


48 


2.969(9) 


0.12 54 


3.401(10) 
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4.461(6) 
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32 


4.488(7) 
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40 


4.546(8) 
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48 


4.543(10) 
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64 


4.579(14) 



0.09 80 5.399(13) 
0.08 92 6.485(22) 
0.07 110 7.917(25) 
0.06 140 9.892(35) 
0.05 24 10.194(7) 
0.05 32 11.108(9) 
0.05 40 11.682(11) 
0.05 48 12.060(13) 
0.05 56 12.279(16) 
0.05 64 12.492(17) 
0.05 80 12.668(21) 
0.05 120 12.892(29) 
0.05 160 12.899(42) 
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C. The thermodynamic Casimir force 

We have computed the thermodynamic Casimir force per area and its first and 
second partial derivative with respect to hi for (0, +) boundary conditions for the 
thicknesses Lq = 8.5, 12.5, and 16.5. To this end we have simulated films of the 
thicknesses Lq = 8, 9, 12, 13, 16 and 17. For most of the simulations we have used 
L = 32 for Lq = 8 and 9, L = 48 for Lq = 12 and 13, and L = 64 for Lq = 16 and 
17. The correlation length of the film displays a single maximum at a temperature 
slightly below the critical temperature of the bulk system. The correlation length 
at the maximum is at most by one per mille larger than at the critical point of 
the bulk system. Therefore our choice of L should ensure that finite L effects of 
the energy per area and its first and second partial derivative with respect to hi 
can be safely ignored. At /3-values that are much smaller or larger than Pc we have 
used smaller values of L. Throughout we have checked that L > IOC, fum is fulfilled 
with a clear safety margin. For Lq = 8 and 9 we have simulated at 85 values of 
the inverse temperature in the range 0.25 < /3 < 0.5, for Lq = 12 and 13 at 124 
values in the rage 0.3 < /3 < 0.42, and for Lq = 16 and 17 at 112 values in the rage 
0.34 < P < 0.406. The difference between neighboring /3-values is adapted to the 
problem: It is the smallest close to (3c- We performed 10^ update cycles for Lq = 8, 
9, 12 and 13 and 2 x 10^ update cycles for Lq = 16 and 17 for each value of /3. 
In total these simulations took about 10 years of CPU time on a single core of a 
Quad-Core AMD Opteron(tm) Processor 2378 running at 2.4 GHz. 

Using the estimates of the energy per area obtained from these simulations we 
have computed the thermodynamic Casimir force per area as discussed in section IVl 
In figure m we have plotted — Lp g^jA/ex as a function of t{LQ^eff/ioY^'^-i where we 
have used -Z^o,e// = Lq + Ls with Lg = 1.43 obtained above in section IVI A[ We do 
not show statistical errors in figure HI since they are comparable with the thickness 
of the lines. The curves for Lq = 8.5, 12.5 and 16.5 fall quite nicely on top of each 
other. Only for x ^ —7, in the low temperature phase, we see a small discrepancy 
between the result for Lq = 8.5, 12.5 and 16.5, which might be attributed to analytic 
corrections. We conclude that we have obtained a good approximation of the finite 
size scaling function 9{q^+). 

Throughout ^(o,+) is positive, which means that the thermodynamic Casimir force 
is repulsive. The scaling function ^(o,+) has a single maximum. We have determined 
the position of this maximum from the zero of AE. We find Pmax = 0.39069(2), 
0.389443(10), and 0.388874(6) for Lq = 8.5, 12.5, and 16.5, respectively. It follows 
xt,max = tmaxiLQ^eff/^o)'^" = "1.184(13), -1.175(11), and -1.174(10) for Lq = 8.5, 
12.5, and 16.5, respectively. The error bar includes the uncertainties of Pmax, Lg and 
u. Note that the results obtained from the three different thicknesses are consistent. 
Next we have determined the value of the scaling function at the maximum. We 
get -Ll^fjAfex{xt,^ax) = 0.567(4), 0.566(3), and 0.564(3) for Lq = 8.5, 12.5 and 
16.5, respectively. The error is dominated by the uncertainty of Lg. The results 
obtained from the three different thicknesses are consistent. As the final result we 
take the one obtained from Lq = 16.5: 

Xt,max = -1-17A{10) , ^(0,+),„^ax = 0.564(3) . (74) 

At the critical point of the bulk system, the finite size scaling function assumes 



23 



-20 -10 10 20 

t(Lo,eff/^o) 



FIG. 4. We plot —Lq ^jj^Afex as a function of i(i^o,e///Co)^^'^ for (0, -|-) boundary 
conditions for the thicknesses Lq = 8.5, 12.5 and 16.5. To this end, we have used 
Lo,e// = Lo + Ls with = 1.43, Co = 0.2282 and = 0.63002. 



the value 

e(o,+)(0) = 0.497(3) (75) 

where the error is dominated by the uncertainty of Lg. This results can be compared 
with 6'(o,+)(0) = 0.33, 0.416 and 0.375(14) obtained by using the e-expansion, and 



Monte Carlo simulations of the Ising model j6l|. Similar to the case of (-I-, +) and 
(-I-, — ) boundary conditions [30], we see a large deviations of the results of Krech 
from ours. 

In figure |5] we compare the finite size scaling function of the thermodynamic 
Casimir force per area for (0, -|-) boundary conditions with those of (-I-, -|-) and 



(-I-, — ) boundary conditions that we have obtained in [30 . 

In the high temperature phase and around the bulk critical point, the absolute 
value of ^(0,-1-) is smaller than that of while in the low temperature phase for 

Xf ~ —1.1 it becomes larger. The value of ^^(o,+) is much smaller than that of _) 
throughout. 



As discussed in ref. [6l(|, see in particular eq. (3.6) and the Appendix A of ref. 
6l| . in the mean- field approximation there is a simple relation between the scaling 
functions &{+,-) and ^(+,o)- For {+, — ) boundary conditions, the magnetization van- 
ishes in the middle of the film. Hence, ignoring fluctuations, a film of the thickness 
2Lq with (-I-, — ) boundary conditions is composed of two films of the thickness Lq, 
where one has (-I-, 0) and the other (0, — ) boundary conditions. Furthermore (0, -|-), 
(-I-, 0) and (0, — ) boundary conditions are equivalent. Therefore 

OMF,(o,+){xt) = 2-''eMF,(+,-){'2}/^Xt) . (76) 
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FIG. 5. We plot our result for the finite size scaling function ^(o,+) along with those for 



\+^+) and ^'(+,-) obtained in ref. 



For less than four dimensions one expects deviations from this relation. Indeed for 
the Ising bulk universality class the ratio of Casimir amplitudes 

= 16(l-0.481e + ...) (77) 

^(0,+) 

obtained by using the e-expansion Glf clearly differs from 2^"*^ = 16(1 — 0.6931. ..e + 
...) obtained from eq. f l76|) . Note that the Casimir amplitude is given by 2A(b^^62) = 
6{bi,b2){^)- For two dimensions one obtains from conformal field theory j62| 



A(+.-) _ 23 
A(+,o) 2 



(7J 



which is almost 3 times as large as the factor 4 predicted by eq. (I75]l . 

Taking our numerical data, we find for > 0, this means in the high temperature 
phase, 9(^+fi){xt) ~ 0.7 x 2^^6'(_|_^_)(2-^/°-^^°°^xt), while in the low temperature phase, 
one gets ^(+,o)(a;t) ^ 2-^e(^+^_){2^/^-'^^^^^Xt) - 0.3 in the range -10 < Xt < -3. This 
means that eq. ( 176|) does not provide a quantitatively accurate relation between the 
scaling functions 9(^+fi){xt) and ^^(+ _)(a;t) in the three dimensional case. 

The most striking observation is that in the high temperature phase d{o,+) decays, 
with increasing xt, much faster to zero than 6{+,+) and ^(+,-) do. This behavior can 
be explained by using the transfer matrix formalism. For a discussion of the transfer 
matrix formalism applied to the problem of the thermodynamic Casimir effect see 
section IV of [30|]. In terms of eigenvalues and eigenvectors \a) of the transfer 



25 



matrix the thermodynamic Casimir force per area can be written as 



T Casimir j- n v— \ / i\ \ \ \ \ V'^/ 





h 


a) (62 


a) 


Eaexp(-m«0 (61 


a 


){b2 


a 


) 



where = '^a = — ln(Aa/Ao). Note that here m is a mass and should not be 
confused with the magnetization. We assume that the eigenvalues are ordered such 
that Xa > A/3 for a < (3, where a, (3 are positive integers or zero. The states 
and 1 62) are defined by the boundary conditions that are applied and I = Lq + 1. 
For 3> the right side of eq. fl7^ is dominated by the contribution from the 
state |1) and therefore 

^(fei,fe2)M) ~ -mH^ exp{-ml)C(hi)C{h2) (80) 
where we have identified 1/^ = m = mi and have defined 

The state |0) is symmetric under the global transformation Sx — for all x in a 
slice. Instead, |1) is anti-symmetric and therefore C = C{+) = — C(— ). It follows 

^(+,+)(m/) = -%,^)(m/) = -C^ mH^ exp{-ml) (82) 

for sufficiently large values of ml. Since Xt = t[l/^oY^'^ — {niiy^'^ it follows 

^^(+,+)(xt) = -9^+,-)ixt) = exp(-xn (83) 

for sufficiently large values of Xt. In the case of free boundary conditions the bound- 
ary state |0) is symmetric under the global transformation Sx — )■ —Sx- Therefore 
vanishes and therefore 

C(0) = . (84) 

Next we have studied the first derivative of the scaling function with respect 
to hi. In figure E] we have plotted — Lq ^ (-^o,e/ / / lex,nor,o) ^olr ^ function of 
t{LQ f,ff/^QY^'^. We do not give error bars, since the statistical error is of similar size 
as the thickness of the lines. We find that the data for Lq = 8.5, 12.5 and 16.5 fall 
quite nicely on top of each other. The small discrepancies that are visible for large 
absolute values of Xt might be attributed to analytic corrections. We conclude that 
our numerical results provide a good approximation of the finite size scaling function 

9'ixt) = _ w/"g read off from figure [6] that 6' is negative throughout 

and has a single minimum. 

We have determined the location of this minimum by searching for the zero of 
We find Pmin = 0.38403(3), 0.38577(2) and 0.38645(2) for Lq = 8.5, 12.5 and 
16.5, respectively. This corresponds to Xt^min = 1-473(18), 1.333(18), and 1.296(24). 
Here we have taken into account the errors of Pmin, Lg and u. In particular for 
Lq = 16.5 the error of Pmin clearly dominates. The results for Lq = 12.5 and 16.5 are 
consistent. As value of the derivative of the scaling function we obtain —0.697(13), 
-0.696(14) and -0.688(13) for Lq = 8.5, 12.5 and 16.5, respectively. Note that in 
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FIG. 6. We plot y = -Ll^fj{Lo^eff/lex,nor,o)-^'^^^§if^ as a function of tiLo^eff/Co)^^" 
for (0,+) boundary conditions for the thicknesses Lq = 8.5, 12.5 and 16.5. To this end, 
we have used Lo^eff = Lq + Lg with Lg = 1.43, = 0.2282 and = 0.63002. 

all cases about half of the error is due to the uncertainty in lex,nor,o = 0.213(3). The 
results for the different lattice sizes are consistent within the quoted errors. We 
conclude 

= 1.30(5) , Cn = -0.69(2) . (85) 

Assuming that C{hi) is an analytic function and the finite size scaling behav- 
ior ( 135|) of the thermodynamic Casimir force per area we arrive at 

e'ixt) = Bxf-^' exp(-a;n (86) 

for Xt ^ 0. Matching our numerical data for Lq = 16.5 at Xt ~ 10 with eq. ( 186|) we 
arrive at B = —0.85(5), where the error is estimated by comparing with the result 
obtained from Lq = 12.5. 

Next we have studied the second derivative of the scaling function with respect 
to hi. To this end, in figure [7] we have plotted -L|j ^ ^ ^ (-^o,e///^ex,nor,o) "^^''i ^ ^Jr 

as a function of t(Lo,e///^o)^'''^- For Lq = 16.5 we have plotted the statistical error, 
which we have not done for Lq = 8.5, 12.5 to keep the figure readable. Within 
our statistical accuracy, the curves for the three different thicknesses fall on top of 
each other. It seems that 6" is positive for all values of the scaling function. Likely 
the negative values found for large \xt\ and Lq = 16.5 are just an artifact due to 
statistical fluctuations. The function displays a single maximum that is located at 

Xt,rmn = "1.9(2) , = "0.39(2) . (87) 
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t(Lo,eff/^o) 

FIG. 7. We plot y = -L^ ^^^(Lo,e///«ex,nor,o)"'^'^i as a function of t(Lo,e///eo)'/" 

for (0,-|-) boundary conditions for the thicknesses Lq = 8.5, 12.5 and 16.5. To this end, 
we have used -Z^o.e// = Lq + Ls with Ls = 1.43, = 0.2282 and v = 0.63002. 



In figure E] we have plotted 6'(o,+), ^(o+) ^'^d O'l^j^y To this end we have used the 
results obtained for Lq = 16.5. We find that the shape of 6'(o+) is quite similar to 
that of 6'(o,+). In particular, for xt — )■ oo, both 6'(o,+) and O'^qj^-^ approach zero much 
faster than 9'^^ Therefore already for an infinitesimally small positive value of , 
the crossover scaling function X/^J, taken as a function of Xt-, has a minimum 
in the high temperature phase. 

In order to check the range of applicability of the Taylor-expansion, and to study 
the crossover beyond the Taylor-expansion, we have simulated films with {hi,+) 
boundary conditions and the thicknesses Lq = 8 and 9 at the values hi = 0.03, 
0.06, 0.1 and 0.2 of the external field at the boundary. Our results along with that 
for -|-) corresponding to hi = jS obtained in |30[ are plotted in figure O For 
hi = 0.03 there is a minimum of the thermodynamic Casimir force per area in 
the high temperature phase. Its absolute value is about one third of the value of 
the maximum in the low temperature phase. The thermodynamic Casimir force 
changes sign at /3 ~ 0.384, which is slightly smaller than Pc- Going to larger values 
of hi the position of the minimum changes only little and the absolute value of the 
minimum increases. On the other hand, the value of the maximum is decreasing 
with increasing hi. For hi = 0.2, the maximum has vanished. 

The authors of 15|] show in figure 9 of their paper Monte Carlo data obtained 
by O. Vasilyev jl6|for the three-dimensional Ising model and the film thickness 
Lq = 10. There is nice qualitative agreement with our results given in figure |9l 

We have compared the results for the thermodynamic Casimir force per area 
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FIG. 8. We plot 0(0,+)) ^(o +) ^(o +) ^ ^ function of xt- 

obtained by simulating at hi = 0.03, 0.06, 0.1 and 0.2 for Lq = 8.5 with those 
obtained by the Taylor-expansion around /ii = up to second order in hi. We find 
that for hi = 0.03 the results almost agree within the statistical error. Still for hi = 
0.06 the Taylor-expansion to second order resembles the true result quite well. The 
largest discrepancy is found for the value of the maximum of the thermodynamic 
Casimir force per area. It is overestimated by about a factor of 1.24. As one 
might expect, the result of the Taylor-expansion becomes increasingly worse with 
increasing hi. In particular it does not reproduce that for large values of hi the 
maximum of the thermodynamic Casimir force per area disappears. 

Given our results for various thicknesses Lq at hi = 0, we conclude that the 
results for Lq = 8.5 provide already a quite good approximation of the scaling 
limit. In particular we are confident that the qualitative features of the crossover 
discussed here still hold in the scaling limit. In particular we conclude that for 
Xh, ^ 0.03[(8.5 + 1.43)/0.213]°-^249 ^ q ^ scaling function e{xt,Xh,) is still well 
described by the Taylor-expansion around Xhi = to second order. 

From figure M we can read off that the thermodynamic Casimir force can also 
change sign as a function of the thickness Lq for fixed values of hi and the temper- 
ature. In general both Xt = t[-Lo/^o]^* and Xf^ = hi[LQ/lex,nor,o]^''^ depend on the 
thickness Lq. Therefore, for simplicity let us consider the bulk critical temperature, 
where Xt = for any thickness of the film. For small Lq the scaling variable Xh^ is 
small and therefore the thermodynamic Casimir force is close to the case Xh^ = 
and is therefore repulsive. As Lq increases, Xhi increases and therefore Q{0,XhJ 
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FIG. 9. We plot — A/e^. for boundary conditions as a function of the reduced 

temperature j3c — (3- The thickness of the film is Lq = 8.5 throughout. 



decreases. We read off from figure [9] that 9(0,a;/iJ ~ for x/^ ~ 1. With further 
increasing Lq, the thermodynamic Casimir force becomes attractive. 



1. Approach to the hi ^ oo limit 



For sufficiently large values of Xh-^ = hi^L^ j lex,noTfiY'^^ we expect that corrections 
to the — )■ oo limit can be described by replacing Lq by L^^^ff = Lq + Ls, where 



In figure [TO] we have plotted our results for hi = 0.2 and Lq = 8.5 and Lq = 16.5. 
First we use L, = 1.9 that we had obtained in ref. for (+, +) boundary 



conditions and second Lg = 1.9 + 1.294, where we have added Al ex, nor {0.2, (3 c) 
obtained in section IVIBI above. For comparison we give the result obtained for 
Lq = 16.5 and (+, +) boundary conditions, using Ls = 1.9. In the case of Lq = 8.5 
the matching with the (+, +) result is somewhat improved by using Ls = 1.9 + 1.294 
instead of Ls = 1.9. While the value of the minimum is clearly improved, the 
matching of the curve with that for (+, +) boundary conditions deep in the high 
temperature phase is not. In contrast, for Lq = 16.5, using Ls = 1.9 + 1.294 instead 
of Ls = 1.9 clearly improves the matching of the curve for hi = 0.2 with that for 
hi = P in the whole range of Xt that is considered. 

We conclude that for Lq ^ 10 Alex, nor {hi, (3 c) using Ls = 1.9 + Alex, nor {hi, /3 c) 

clearly improves the matching with the (+, +) scaling function. It would be desir- 
able to check this by simulations for smaller values of hi. However this would be 
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FIG. 10. We plot Ll^^j^jAfex as a function of tiLo^effHoY'" for hi = 0.2 and Lq = 8.5 
and 16.5 using Lg = 1.9 or Ls = 1.9 + 1.294. For comparison we give the corresponding 
curve for Lq = 16.5 and (+, +) boundary conditions using Ls = 1.9. For a discussion see 
the text. 

quite expensive, since already for hi = 0.15 we would need to simulate a thickness 
Lo ^ 30. 

VII. SUMMARY AND CONCLUSIONS 

We have studied the crossover behaviors of a surface of a system in three- 
dimensional Ising universality class from the ordinary to the normal or extraor- 
dinary surface universality class. To this end, we have simulated the improved 
Blume-Capel model on the simple cubic lattice. In particular we have studied films 
with various boundary conditions applied. Improved means that corrections to fi- 
nite size scaling oc L^'^ have a vanishing amplitude, where Lq is the thickness of 
the film and u = 0.832(6) [36j is the exponent of leading corrections. This property 
is very useful in the study of films, since corrections oc Lq ^ due to the surfaces are 
expected js^] and fitting data it is difficult to disentangle corrections with similar 
exponents such as u and one. Mostly we have simulated films with (0, +) bound- 
ary conditions. This means that at one surface we apply free boundary conditions, 
while at the other surface the spins are fixed to Studying the magnetization of 
the slice at the surface with free boundary conditions, at the bulk critical point, of 
films of a thickness up to Lq = 64 we arrive at the estimate i/hi = 0.7249(6) for the 
renormalization group exponent of the external field at the surface for the ordinary 
surface universality class. This estimate is at least by a factor of 5 more accurate 
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than those previously given in the hterature. The authors of [51| quote an error 
that is only 2.5 times larger than ours, however the deviation between our and their 
estimate is about 6 times larger than the combined errors. For details see table 
[H We have studied the magnetization profile in the neighborhood of the surfaces 
for both the ordinary as well as the normal surface universality class. The data 
are consistent with the theoretically predicted power law behavior. This study also 
allowed us to determine the extrapolation length /ex for free boundary conditions as 
well as symmetry breaking boundary conditions for various values of the external 
field hi at the surface. Corrections to scaling oc Lq^, which are due to the surfaces 
of the film can be expressed by an effective thickness -Z^o,e// = -^^o + Lg, where Lg 
depends on the details of the model. Our numerical results confirm the hypothesis 
that Lg = lex,i + hx,2, where lex,i and lex,2 are the extrapolation lengths at the two 
surfaces of the film. 

Next we have studied the thermodynamic Casimir force in the neighborhood of 
the bulk critical point in the range of temperatures where it does not vanish at the 
level of our accuracy. First we have simulated films with (0, +) boundary conditions 
and the thicknesses Lq = 8.5, 12.5 and 16.5. Taking into account corrections by 
replacing Lq by L^^eff-, the behavior of the thermodynamic Casimir force and its 
first and second derivative with respect to hi follows quite nicely the predictions 
of finite size scaling. Hence our data allow us to compute good estimates of the 
finite size scaling functions 6^(0,+), ^(o+)) and 6'(o_,_)- Next we have computed the 
thermodynamic Casimir force per area for the thickness Lq = 8.5 at the finite values 
hi = 0.03, 0.06, 0.1 and 0.2 of the external field at the boundary. We find that the 
Taylor-expansion of the thermodynamic Casimir force up to the second order in hi 
around hi = still describes the full function well at hi = 0.03 which corresponds 
to the value = /ii[i^o/^o,ea:,nor]^''i ~ 0.5 of the scaling variable of the external 
field at the boundary. Finally, we have studied the approach of the thermodynamic 
Casimir force to the limit /ii — )■ oo. We find that by using I/o,e// = Lo + Ls{hi), the 
corrections to this limit are well described for Lq ^ 10[Ls{hi) — Ls{/3c)]- 

Based on exact results for stripes of the two-dimensional Ising model [13], mean- 



field calculations [15| and preliminary Monte Carlo results for the Ising model [16 



on the simple cubic lattice one expects that for certain combinations of the external 
fields hi, h2 and the thickness of the lattice Lq, the thermodynamic Casimir force 
changes sign as a function of the temperature. Also for certain choices of the external 
fields hi, h2 and the temperature, the thermodynamic Casimir force changes sign as 
a function of the thickness of the film. Here we confirm these qualitative findings. 
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